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1  Introduction 


Mature  and  accurate  computational  models  exist  for  many  oceanographic  applications;  how¬ 
ever,  the  required  model  inputs,  in  particular  boundary  condition  data,  are  not  directly 
available  in  many  cases  of  practical  interest.  The  data  that  does  exist  tends  to  be  at  the 
wrong  locations,  interior  to  the  spatial  domain  of  interest  and  not  on  the  boundary.  This 
incompatibility  between  the  required  model  inputs  and  the  available  data  represents  a  sig¬ 
nificant  obstacle  to  using  observational  data  to  improve  model  predictions. 

In  this  study  the  use  of  observational  data  in  model  calculations  of  the  wave  spectrum  in 
the  near-shore  region  is  examined.  The  model  used  is  the  SWAN  model  [18,  4],  a  near-shore 
wave-spectrum  model,  which  implements  conservation  of  the  wave-action  spectral  density  in 
shallow  water.  The  model  can  include  the  effects  of  variable  bathymetry,  currents,  bottom 
friction,  depth-induced  breaking,  wind  input  and  whitecapping.  The  problem  to  be  examined 
is  the  use  of  observations  at  interior  points  in  the  region  of  interest  to  estimate  the  appropriate 
incident- wave  spectrum  boundary  conditions  and/or  the  bottom  friction  coefficient  under 
steady-state  conditions. 

The  observational  data  to  be  examined  will  include  satellite-based  SAR  imagery,  as  well 
as  in-situ  observations  of  the  wave  directional  spectrum.  For  SAR  data  assimilation,  a  fully 
nonlinear  mapping  approach  [10,  12]  is  used  to  estimate  the  SAR-image  spectrum  from  the 
wave  spectrum.  In  the  mapping  from  the  wave  spectrum  to  the  SAR  image,  some  of  the 
wave  information  in  lost;  energy  at  high  wavenumbers  in  general,  and  at  high  azimuth  wave 
numbers,  in  particular,  is  lost.  This  effect  will  be  apparent  in  some  of  the  results. 

The  variational  approach  adopted  here  follows  that  of  Bennett  [2]  in  its  derivation,  and 
that  of  LeDimet  &  Talagrand  [14]  in  application.  An  objective  function  is  defined  which 
quantifies  the  error  in  the  model  prediction  of  the  observation  data.  To  ensure  the  uniqueness 
of  the  result,  the  cost  function  includes  a  term  which  ‘penalizes  the  estimated  model  inputs 
[3[.  An  adjoint  SWAN  model  equation  is  derived  which  has  as  an  input  the  difference  between 
the  model  prediction  and  the  observation  data.  From  the  adjoint  solution  the  gradient  of  the 
objective  function  with  respect  to  the  model  inputs,  the  incident  wave  spectrum  and/or  the 
bottom  friction  coefficient  are  calculated.  The  gradient  is  used  together  with  a  conjugate- 
gradient  minimization  algorithm  to  determine  the  model  inputs  which  minimize  the  objective 
function. 

In  the  next  section,  the  SWAN  model  and  the  SAR-image  spectrum  model  is  described. 
Then  the  objective  function  is  defined  and  the  adjoint  SWAN  model  is  derived,  along  with 
expressions  for  the  gradient  of  the  cost  function  with  respect  to  the  model  inputs.  After  that, 
the  algorithm  is  applied  to  synthetic  data  and  real-world  data:  both  in-situ  spectrum  mea¬ 
surements  and  SAR-image  data.  In  the  real-world  data  cases,  the  results  for  the  estimated 
wave  field  are  evaluated  by  comparison  to  independent  data. 

2  The  Models 

The  assimilation  procedure  relies  on  two  models,  one  which  governs  the  wave  spectrum  and 
another  which  relates  the  wave  spectrum  to  the  spectrum  of  the  SAR  image. 
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2.1  The  Wave  Spectrum  Model 

The  SWAN  model  [18,  4]  is  a  near-shore  wave-action-balance  model  which  can  predict  the 
evolution  of  the  wave  spectrum  in  coastal  regions.  The  wave-action  spectral  balance  is 
expressed  as 


dN_ 

~dt 


t  V  ■  (C N) 


S{N) 

a 


(i) 


where  x  and  y  are  spatial  position  variables,  and  a  and  9  are  wave  frequency  and  direction 
variables  for  the  action  spectrum,  V  =  (^,  J|),  and  C  =  (Cx,Cy,Ca,Ce)  represents 

the  wave-energy  propagation  velocities  in  physical  and  spectral  space.  Below,  the  vectors 
x  =  ( x,y )  and  s  =  (<j,  0)  will  be  used  to  represent  spatial  and  spectral  position,  respectively. 
iV(x,  s,  t)  is  the  action  spectral  density  defined  as 


N(x,s,t)  —  E(x,s,t)/a, 


(2) 


where  E  is  the  energy  spectral  density  and 

a(k)  =  \J gk  tanh  kh.  (3) 

is  the  intrinsic  radian  frequency,  where  g  is  the  gravitational  acceleration,  h(x)  is  the  water 
depth,  and  k  -  2ir/\  is  the  wavenumber  (A  is  the  wave  length).  The  x  and  y  -direction 
components  of  the  wave-propagation  velocities  are  given  by 

Cx  =  U  +  Cgcosd,  Cy  =  V  +  CgsmO,  (4) 

where  U(x,t )  and  V(x,t)  are  the  x  and  y  direction  current  velocities,  specified  as  inputs  to 
the  problem,  and  the  wave  group  velocity  is 

c’  =  Ui  tanh“(1  +  ^*)  '  <5) 

The  other  two  velocities  Ca  and  Cg  are  energy  propagation  velocities  in  the  spectral  domain 
caused  by  depth  and  current  variations.  They  are  most  easily  defined  in  terms  of  the  apparent 
frequency  fi  (as  seen  by  a  stationary  observer),  which  includes  a  doppler  shift  induced  by 
the  current 


Q  =  a  -f  k  (U  cos  d  +  V  sin  9) . 


(6) 


The  spectral  propagation  velocities  are  given  by 


Ca  = 


dQ 
dt  ’ 


'dQ  .  .  dtl 
— —  sin  6  +  —  cos  9 
ox  oy  . 


(7) 


The  source  term  on  the  right-hand  side  of  (1)  is  described  in  detail  in  [18]  and  includes 
the  effects  of  wind  growth  and  energy  transfer  in  the  spectrum  due  nonlinear  wave-wave 
interactions  (resonant  triad  and  quartet  interactions).  Significant  additional  contributors 
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to  the  source  term  are  various  processes  by  which  wave  energy  is  dissipated.  These  include 
whitecapping,  bottom  friction  and  depth-induced  breaking.  Energy  loss  due  to  whitecapping 
is  dependent  on  wave  steepness  and  is  parameterized  following  the  pulse-based  model  of  [8], 
adapted  by  VVAMDI  Group  (1988), 


SdlW  =  -Yal  E 

k 


(8) 


where  a  is  the  mean  intrinsic  frequency  and  k  is  the  mean  wavenumber.  T  is  a  steepness- 
dependent  coefficient  estimated  by  examining  fully  developed  conditions.  The  general  form 
of  the  bottom  friction  source  term  is 


(9) 


where  Q,  is  a  spatially  varying  dissipation  coefficient  that  can  be  parameterized  as  a  function 
of  the  bottom  roughness,  integral  spectral  parameters,  or  the  frequency  and  direction  of  a 
wave  component  (see  8],  [9],  [5],  [15]  ).  For  depth-induced  breaking,  SWAN  essentially  pre¬ 
serves  the  spectral  shape  (dissipates  the  spectrum  uniformly)  and  uses  the  simple  saturated 
breaker  criterion  with  an  allowed  variability  for  bottom  slope 


(10) 


where  S^br.tot  is  the  rate  of  total  wave  energy  dissipation  and  is  determined  following  Battjes 
and  Janssen  (1978). 

This  set  of  equations  can  be  solved  for  the  action  spectrum  for  spatial  region  77  subject  to 
appropriate  boundary  conditions  on  <977.  For  portions  of  the  wave  spectrum  with  propagation 
velocities  which  carry  energy  into  77,  the  ‘incident’  wave  spectrum  JV(x,  s,  t)  on  the  boundary 
(777  must  be  specified;  for  the  portions  of  the  wave  spectrum  for  which  the  propagation 
velocities  carry  the  energy  outward  from  77,  an  ‘outflow’  condition  is  used  (continuity  of 
the  action  flux  at  <977).  In  the  spectral  domain,  for  most  practical  implementations,  the 
spectral  density  is  required  to  vanish  on  the  upper  and  lower  frequency  (cr)  boundaries;  this 
condition  is  satisfied  by  locating  the  a  boundary  far  from  the  energy-containing  region  of  the 
spectrum.  The  boundary  conditions  in  6  are  that  the  spectrum  is  periodic.  In  addition  to 
these  boundary  and  initial  conditions,  complete  specification  of  the  mathematical  problem 
requires  the  the  bathymetry  fi(x)  and  current  field  U (x,  t ),  V(x,  t)  to  be  prescribed  for  spatial 
region  77. 

2.2  The  SAR-lmage  Spectrum  Model 

The  SAR  -image  spectrum  can  be  related  to  the  wave  spectrum  through  either  a  fully 
nonlinear  mapping,  or  a  less-accurate  but  simpler,  quasi-linear  model.  These  are  described 
in  the  next  sections. 
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Nonlinear  SAR  Model 

A  model  which  relates  the  SAR-image  spectrum  to  the  ocean-wave  spectrum  has  been  derived 
by  Hasselmann  &  Hasselmann  [10]  and  Krogstad  [12].  As  summarized  by  Lyzenga  [13],  the 
SAR-image  spectrum  is  given  by 

/(k)  =  i//G(r'W""kr*  (H) 

where  r  is  a  position  vector  in  correlation  space.  ka  is  the  azimuth  (along-track)  wavenumber 
component  of  k,  the  wave  vector.  The  term  in  the  integrand  is  defined  as 

G{ r,  ka)  =  [prr(r)  +  /( r,  ka)f(-r,  -fca)]  e~^C(r),  (12) 


where 


—  f  ika  [Prv(b)  PrtXOj  j 
C(r)  =  Pw(0)  -  Pw(r), 


(13) 

(14) 


and  the  spatial  correlation  functions  prr  and  pvv,  as  well  as  the  cross-correlation  function 
prv,  are  defined  as 


Prr (r) 

=  /  j  \Tr(k)\2  E(k)eikrdk, 

(15) 

Pw(  r) 

=  j  j  |Ti,(k)|2  E(k)elk'Tdk,  and 

(16) 

Prv{  r) 

=  f  jTr(k)T;<k)E(k)e*rdk. 

(17) 

Here, 


Fi 

T,(k)  =  — —  (gkT  sin  $,  —  ia2  cos  Rj  , 


(18) 


where  R  is  the  range  from  the  sensor  to  the  scene  center,  V  is  the  sensor  velocity,  a  is  the 
wave  frequency  from  (3),  Ay  is  the  range  (cross-track)  component  of  k  and  0i  is  the  incidence 
angle.  The  other  function  is  given  by 

Tr(k)  =  k  [mft(k)  +  imt( k)j ,  (19) 

where  k  =  |k|.  The  ‘hydrodynamic’  modulation  of  the  radar  reflectivity  is  given  by 

m‘(k,  =  5^M\^cos2(s-^'  (20) 

where  0r  is  the  radar  look  direction  and  6  is  the  wave  propagation  direction;  m. o  =  7.5  for 
vertical  polarization  and  mo  =  12.5  for  horizontal  polarization.  The  ‘tilt’  modulation  is  given 

by 


mt(k)  =  (5  cot  6i  +  4 tan  0,:  —  4  7P  sin  0,)  cos(0  —  9r), 


(21) 
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with 


lv 


1 

cos  6i  +  n 

n  cos  $i 

n  cos  $i  +  1  1  T  sin2  0* 


for  horizontal  polarization, 
for  vertical  polarization, 


and  for  sea  water  n  «  9. 


(22) 


Quasi-Linear  SAR  Model 

As  is  evidenced  by  (11)  through  (22),  the  full  nonlinear  mapping  between  the  ocean-wave 
spectrum  E  and  the  SAR-image  spectrum  /  is  quite  complicated.  A  simpler,  quasi-linear 
approximation  [10]  can  be  expressed  as 

7(k)  »  Mg( k)  E{ k)  +  M's(- k)  E(- k)  (23) 

where 

M's (k)  =  i  |Tr(k)  +  ika%( k)|2 e-jk"'T"  (24) 

where  Tr( k)  and  Tv(k)  are  as  specified  above,  and  cr2  =  pw(0).  For  use  in  the  assimilation 
procedure,  this  can  be  re-written  in  polar  s  =  (a,  6)  coordinates  as 

f(s)  =  Ms(s)  E(s)  +  M${s~)  E(s~),  (25) 


where  s  =  (a,  9  —  tt)  and 


Ms(s)  =  k—M's(k).  (26) 

For  the  results  shown  herein,  the  full  Hasselrnann  &  Hasselmann  model  (11-22)  was  used 
for  forward  prediction  of  the  SAR-image  spectrum,  while  the  adjoint  SAR  model  used  in  the 
inversion  procedure  was  based  on  the  quasi-linear  model  shown  in  (24). 


3  Assimilation  Methodology 

In  this  section  the  basis  for  the  variational  data  assimilation  approach  is  developed.  The 
problem  setup  will  follow  the  ‘weak  constraint’  formulation  of  Bennett  [2],  but  the  strong 
constraint  limit  will  be  taken,  yielding  an  approach  similar  to  that  of  Le  Dimet  &  Talagrand 
[14].  The  model  inputs  being  estimated  are  penalized  in  the  objective  function  to  ensure 
uniqueness  [3],  An  objective  function  is  defined  which  is  a  positive-definite  measure  of  the 
difference  between  a  set  of  observations  and  the  model  predictions,  as  well  as  the  fit  of  the 
models  to  the  data.  This  objective  function  is  to  be  minimized  by  adjusting  the  S WAN- 
model  inputs.  For  the  purposes  of  this  study,  only  stationary  conditions  will  be  considered, 
and  the  model  inputs  are  taken  to  be  the  boundary  conditions  (incident-wave  spectra)  and 
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the  bottom  friction  coefficient. 

In  what  follows,  we  first  define  the  objective  function  to  be  minimized  (J).  We  then 
determine  the  conditions  for  a  minimum  in  J,  the  Euler-Lagrange  equations.  This  is  fol¬ 
lowed  by  the  relations,  deriver  from  the  Euler-Lagrange  equations  used  in  the  assimilation 
algorithm.  The  section  closes  with  the  steps  implemented  in  the  assimilation  algorithm. 


3.1  The  Objective  Function 

Our  goal  is  to  minimize  the  error  in  our  prediction  of  some  observable  quantity  D(x,  s) 
compared  to  a  set  of  observations.  The  observation  data  D  are  taken  to  be  related  to  the 
wave  energy  spectrum  E  by 


D(x,  s)  =  M(s)E(x,  s)  +  M(s-)E(x,s~),=  ME  +  {ME)~, 


(27) 


where  M (x,  s)  is  in  this  case  a  transfer  function.  For  direct  observations  of  the  wave  spec¬ 
trum,  M  =  1  and  ( ME)~  =  0,  while  for  SAR  observations,  M  =  Ms  from  (24),  above.  For 
a  set  of  Nd  observations  at  spatial  locations  x,,  the  error  variance  between  the  predictions 
and  observations  can  be  expressed  as 


1 

N~d 


Di)25Xi  ds  dx, 


where  Dt  is  an  observation  of  D  at  location  xu  5Xi  =  3(x  -  x,)  is  a  function  sampling  D  at 
x,  for  comparison  to  the  data,  and  the  integral  is  over  the  entire  spectrum. 

To  obtain  a  prediction  from  our  coupled  wave-spectrum/SAR -spectrum  model  which 
matches  a  set  of  observations,  we  first  define  an  ob  jective  function  to  be  minimized,  which 
includes  both  the  data  and  the  models  used  to  predict  the  data,  along  with  any  additional 
constraints  to  be  placed  on  the  solution.  In  this  approach,  we  wish  to  minimize  the  error  in 
the  predictions  of  the  observed  data,  as  well  as  the  misfit  between  the  models  and  the  actual 
physics.  The  objective  function  J  is 


J 


+ 


+ 


1  r  r  Nd  _ 

lUs£D- 


Di)26Xt  ds  dx 


1  PC  \  2 

-2Us^\Z[D-  ME -(ME)-}  6, 


Ar 


2  Jb  Js  Lb 


wj-p-  E2  ds  dxb 


V  ■  (C N)  +  Cbf 
1 


E 


-,2' 


ds  dx 


try--—  (Cb  f)2  ds  dx  , 


2  JnJs  Ato 


(28) 


where  the  first  line  is  recognized  as  proportional  to  the  error  variance  between  the  model  and 
the  data,  the  second  and  third  contain  the  mean-square  misfit  of  the  models  for  the  SAR 
spectrum  and  the  wave  spectrum.  The  last  line  contains  positive-definite  measures  of  the 
model  inputs  to  be  estimated,  the  boundary  spectrum  E(xb,  s),  and  the  spatially  varying 
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bottom  friction  coefficient  Cb]  these  are  included  to  ensure  a  unique  solution  [3].  In  (28), 
Ar  is  the  area  of  the  spatial  domain,  (f>m  is  a  constant  with  the  same  units  as  M'2,  Lb  is  the 
length  of  the  boundary,  xb  is  the  position  along  the  boundary,  and  4>j  is  a  constant  with 
the  same  units  as  D.  Here  w  is  a  weighting  factor  controlling  the  contribution  of  the  model 
misfit  to  the  objective  function,  relative  to  the  data,  and  Wb  and  wj  are  weighting  factors 
controlling  the  contributions  of  the  estimated  model  inputs. 

3.2  The  Euler-Lagrange  Equations 

The  objective  function  J  depends  on  D,  E,  Cb  and  the  models  that  relate  them.  The 
best  lit  of  the  model  results  to  the  data  will  be  for  the  boundary  spectrum  and  bottom 
friction  coefficient  which  makes  J  stationary,  i.e.  when  the  first  variation  of  J  with  respect 
to  these  variables  vanishes.  The  first  variation  of  (28)  with  respect  to  D,  E  and  Cb  has  three 
components 


8J  =  8Jd  +  8JE  +  8JCb-  (29) 

The  relations  which  define  the  minimum,  which  are  satisfied  when  8  Jo-,  8JE,  and  5Jcb  vanish 
are  the  Euler-Lagrange  equations. 

The  first  component  of  8J  is  given  by 

/.  Nd  _ 

SJD  =  /  E  \(D  “  A)  +  w[D  ~  ME  ~  (ME)-)}  5Xi  8D  ds  dx,  (30) 

J  TZ  S  j  j 

Since  the  region  over  which  the  integration  occurs  is  essentially  arbitrary,  and  5D  is  in  general 
non-zero,  8 Jo  will  vanish  when 


Nd 


\D  -  ME  - 


Nd 


■>*.  =  -£ -(£>- a)  * 


i=l  i  \ 

The  second  component  of  §J  is,  after  integration  by  parts. 


(31) 


-  ME  -  {ME)~  M8Xi  +  C  •  VT  -  CbfA 


8E  ds  dx 


+ 


j  J  \  wb~p-  E  +  AC  ■  nj  8E  ds  dxb, 


(32) 


where  n  is  a  unit  vector  normal  to  the  domain  boundary,  and 


A  = 


m>m 
Ara  . 


V-  C N)+Cbf 


E 


a 


(33) 


Both  integrals  in  (32)  will  vanish  independently  at  the  minimum  for  J.  Requiring  the  first, 
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area,  integral  to  be  zero  yields 


~  ~  Nd  r 

C  ■  VA  =  CbfA  -wJ2[D-  ME  -  (ME)”]  M5Xl. 

i—  1 


(34) 


and  (32)  becomes 


=  U( 


wb-p  E  +  AC  ■  n  \  5E  ds  dxb, 
Eb  J 


(35) 


where  only  the  integral  along  the  boundary  is  retained.  The  final  component  of  5J  is 

5JCb  =  -  /  /  f/AE  -  wf^Cb)  5Cb  ds  dx  .  (36) 

Ju  JS  \  Aro-  ) 

At  the  minimum  in  J  we  must  have  then 


0  =  I  WbjE  E  +  AC  •  nj  5E  ds  dxfc 


(37) 


and 


0  =  - 


K  -JS 


fAE  -  wj^—lcA  SCb  ds  dx  . 


(38) 


The  Euler  Lagrange  equations,  those  which  define  the  minimum  for  J,  are  (31),  (34), 
(37)  and  (38),  along  with  definition  (33). 


3.3  The  Assimilation  Equations 

The  assimilation  procedure  is  intended  to  find  the  model  inputs  Eb  and  Cb  which  results  in 
a  solution  E(x,  s)  that  satisfies  the  Euler-Lagrange  equations.  The  equations  used  in  the 
assimilation  procedure  include  the  governing  equation  for  E(x,  s), 


V  ■ 


E  Ato 
-Cbf-  +  — —  A, 

O  W0m 


(39) 


and  the  ‘adjoint’  equation  governing  A(x,  s)  which  results  from  combining  (31)  and  (34), 


N, 


c  •  VA  =  C6/A+  ID  -  DA  5 


(40) 


i  1 


The  adjoint  equation  includes  the  model-data  error,  and  the  transfer  function  M  from  the 
sensor  model  (27),  on  the  right-hand  side.  To  obtain  expressions  for  the  gradient  of  J  with 
respect  to  the  model  inputs,  we  use  (35),  which  leads  to 

3J  f  ~ 

— -  =  (j)mwbEb  +  /  AC  •  n  dxb,  (41) 

oEb  Jb 
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for  a  boundary  spectrum  which  is  constant  in  space,  E(xt,s)  =  Eb( s),  and  (36),  which  leads 
to 


(42) 


At  the  minimum  for  the  objective  function  J,  these  gradients  will  be  zero,  thereby  satisfying 


the  Euler-Lagarange  equations. 

Equations  (39)  through  (42),  along  with  the  sensor  model  (27),  constitute  the  full  set  of 
equations  for  our  data  assimilation  problem.  This  is  a  weak  constraint  formulation,  where 
the  governing  SWAN  model  and  the  sensor  model  are  not  required  to  be  exactly  satisfied. 
As  a  result,  (39)  and  (40)  are  coupled  through  forcing  terms  on  their  right-hand  sides,  and 
require  a  simultaneous  solution.  These  equations  can  be  uncoupled  by  allowing  w  — *■  oc, 
which  causes  (39)  to  revert  back  to  (33).  This  results  in  the  ‘strong  constraint’  formulation, 
where  the  forward  models  are  taken  as  exact,  and  all  the  error  is  assumed  to  reside  in  the 
data. 

3.4  The  Assimilation  Algorithm 

For  the  results  presented  below,  wave-spectrum  or  SAR-image-spectrum  observations  at 
one  or  more  points  are  used. The  initial  guess  for  the  wave  spectrum  over  the  entire  region 
is  zero,  including  the  incident  wave  spectrum.  The  adjoint  equation  (40)  is  solved  with 
the  observations  as  input,  and  the  adjoint  solution  is  used  to  calculate  the  gradient  (41), 
(42).  The  gradient  is  used  in  a  conjugate-gradient  procedure  1 16]  to  find  the  incident  wave 
spectrum  and  bottom  friction  coefficient  which  results  in  the  minimum  in  the  cost  function 
J.  In  practice,  the  following  sequence  of  steps  is  followed: 

1.  The  adjoint  solution  is  calculated  using  using  the  error  in  the  most  recent  prediction 
as  input  (i.e.  the  second  term  on  the  right-hand  side  of  equation  40). 

2.  From  the  adjoint  solution,  the  gradient  is  determined  using  (41),  (42). 

3.  The  conjugate  gradient  algorithm  is  used  to  calculate  a  new  estimate  of  the  incident 
wave  spectrum  Ef,  and/or  bottom  friction  C\. 

4.  The  SWAN  model  is  run  with  the  new  inputs  and  a  new  wave-spectrum  prediction  for 
the  region  is  generated. 

5.  The  sensor  model  is  run  with  the  new  spectrum  as  input  and  a  new  prediction  of 
the  data  is  generated.  For  SAR-image  data,  the  sensor  model  is  the  Hasselmann  & 
Hasselmann  fully  nonlinear  model  (11),  while  for  in-situ  data  the  sensor  model  is  a 
unit  transfer  function. 

These  steps  are  repeated  until  the  gradient  of  J  goes  to  zero,  which  occurs  when  the  cost 
function  ceases  to  change.  The  final  prediction  from  the  SWAN  model  is  considered  the 
‘best-fit’  to  the  observation  data. 
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4  Application  of  the  Assimilation  Algorithm 

Here  we  apply  the  assimilation  algorithm  three  ways:  First  we  apply  it  to  in-situ  wave- 
spectrum  observations  to  estimate  the  incident  wave  spectrum  which  yields  a  prediction 
matching  the  observations.  This  is  done  first  for  synthetic  observations  (generated  using  the 
SWAN  model),  and  then  for  actual  field  data.  The  algorithm  is  then  applied  to  SAR-image 
data,  again  being  used  to  estimate  the  incident  wave  spectrum  for  first  simulated  SAR- 
image  spectra  and  then  for  actual  ERS  SAR-image  spectra.  Finally,  initial  results  showing 
estimation  of  bottom  friction  from  ERS  SAR  data  is  presented. 

The  assimilation  procedure  will  be  applied  for  a  region  off  the  Atlantic  coast,  centered 
on  the  US  Army  Corps  of  Engineers  Field  Research  Facility  (FRF)  in  Duck,  NC.  The 
bathymetery  for  the  region  is  shown  in  figure  1  and  extends  277.8  km  North-South,  from 
South  of  Cape  Hatteras  to  North  of  the  mouth  of  the  Chesapeake  Bay,  and  134.4  km  East 
West,  so  that  it  extends  offshore  beyond  the  continental  shelf.  For  the  cases  examined  below, 
the  region  shown  in  figure  L  is  discretized  at  995  m  (E-W)  and  2013  m  (N-S)  resulting  in 
a  130  x  139  computational  grid.  The  origin  of  the  coordinate  system,  the  lower,  left-hand 
corner  of  the  region  in  figure  1.  is  located  it  35CN  76° W.  and  it  extends  from  35°  —  36.5°N 
.  lid  7o  -  71  .5  W  The  v ,.i\"  spectrum  has  25  frequencies  logarithmically  spaced  from  0.042- 
ii  ill  Hz.  and  4s  equallv  spaced  directions  covering 360  .  A  spatially  uniform  wave  spectrum 
at  the  Eastern  incident-wave  boundary’  is  considered  and  so  it  is  desirable  that  the  Eastern 
boundary  be  in  uniformly  'deep'  water  (i  e.  that  none  of  the  waves  are  influenced  by  the 
bottom)  It  is  assumed  that  no  waves  enter  at  the  North  and  South  boundaries.  At  the  shore 
boundary,  the  action  spectral  density  vanishes.  Due  to  the  lack  of  waves  on  the  North  and 
South  boundaries,  the  solution  will  not  be  valid  in  triangular  regions  which  grow  Westward 
from  the  Eastern  corners  of  the  domain.  The  region  chosen  will  therefore  allow  a  portion  of 
the  coastline  roughly  100  150  km  in  extent  to  be  accurately  simulated  with  the  specification 
of  a  constant  incident  wave  spectrum. 

In  practice,  the  SWAN  model  (v40.11)  is  used  for  the  forward  predictions  of  the  wave 
spectrum.  A  modified  version  of  the  SWAN  model  (v30.75)  is  used  to  solve  the  adjoint 
equations.  The  calculation  of  the  cost  function  and  the  observation  term  for  the  adjoint 
equations  is  done  using  a  separate  FORTRAN  code.  The  conjugate  gradient  minimization 
algorithm  and  calculation  of  the  incident  wave  spectrum  are  implemented  in  yet  another 
FORTRAN  code.  The  iteration  scheme  outlined  above  is  carried  out  by  executing  the 
various  codes  under  control  of  a  C-shell  script.  For  the  SWAN  model  predictions,  depth- 
induced  wave  breaking  and  bottom  friction  are  included,  but  not  wind  growth,  white-capping 
or  nonlinear  interactions.  For  the  adjoint  SWAN  solutions,  all  source  terms  are  ignored  in 
order  to  simplify  the  implementation  of  the  adjoint. 

4.1  Estimation  of  the  Wave  Field  from  In-Situ  Observations 

Application  to  Simulated  Wave-Spectrum  Data 

In  order  to  assess  the  accuracy  of  the  assimilation  procedure,  it  will  be  applied  initially  to  a 
synthetic  data  set,  generated  by  a  forward  run  of  the  SWAN  model.  Use  of  synthetic  data 
allows  the  procedure  to  be  evaluated  independent  of  the  accuracy  of  the  SWAN  model  itself. 
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x  (km) 


Figure  1 .  Color  plot  of  the  water  depth  fi(x)  in  meters  for  a  region  roughly  centered 
on  the  USACE  FRF  at  Duck,  NC:  'k,  location  of  the  FRF  8  m  array;  0, 
location  of  NDBC  Buoy  44014;  A,  location  of  NDBC  CHLV2  station. 
The  x  and  y  axes  correspond  to  longitude  and  latitude,  respectively,  and 
the  lower  left-hand  corner  corresponds  to  35°N  76° W 


To  get  an  indication  of  the  results  from  the  procedure  for  real-world  data,  it  will  then  be 
used  with  a  data  set  from  the  FRF  8m  array  directional  spectrum  as  input.  The  resulting 
wave  field  estimate  will  be  compared  to  the  input  data,  and  also  to  data  from  two  NDBC 
buoys  located  in  the  region  of  interest  which  were  not  used  in  the  assimilation  procedure. 
The  locations  for  these  data  are  shown  in  figure  1. 

To  generate  the  synthetic  data  set,  the  SWAN  model  was  run  for  the  area  shown  in 
figure  1.  The  incident- wave  spectrum  was  a  Pierson-Moskowitz  spectrum,  where  wind  speed 
at  10  m  above  the  ocean  surface  (C/jo)  and  wind  direction  ($p) are  the  defining  parameters. 
The  Pierson-Moskowitz  spectrum  is  given  by 


E(cr,6)  =  ag2(2ir)  4/  5  exp 


sech2  \p{6  -  0P)} 


(43) 


where  /  is  wave  frequency,  a  =  0.0081  and  fm  =  O.l3g/U\o  [!]•  Here,  the  directional  spread 
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is  modeled  after  [7]  with  P  =  2.  The  example  to  be  shown  will  be  for  a  case  where  the 
wind  direction,  and  hence  the  dominant  wave  direction,  is  dp  =  165°  counterclockwise  from 
East  arid  U\ o  =  12  m/s.  The  resulting  wave  spectrum  is  shown  in  figure  2a).  For  this 
spectrum,  the  significant  wave  height  Hs  =  3.15  m  and  the  peak  frequency  is  fp  =  0.108  Hz 
when  interpolated  onto  the  SWAN  spectral  grid.  This  spectrum  was  used  as  the  input 
wave  spectrum  along  the  right-hand  (Eastern)  boundary  of  the  region  shown  in  figure  1. 
The  SWAN  model  solution  is  shown  in  figure  2b)  and  2c).  Figure  2b)  shows  the  directional 
spectrum  at  the  location  of  the  FRF  8  m  array.  The  dominant  wave  direction  at  this  location 
has  rotated  to  9V  —  183.3°,  roughly  the  cross-shore  direction,  and  the  directional  spread  of  the 
spectrum  has  been  considerably  reduced  relative  to  that  shown  in  figure  2a).  The  significant 
wave  height  has  been  reduced  to  Hs  —  2.63  m  due  to  energy  dissipation  from  bottom  friction 
and  depth-induced  wave  breaking.  The  resulting  significant  wave  height  field  for  the  region 
is  shown  in  figure  2c).  Here,  the  the  reduction  in  Hs  near  the  upper  and  lower  boundaries, 
due  to  the  boundary  conditions  there,  can  be  seen. 

The  spectrum  from  the  8  m  array  location  was  used  as  synthetic  data,  and  assimilation 
into  the  SWAN  model  was  accomplished  using  the  methodology  described  above.  The  initial 
guess  was  assumed  to  be  a  zero  incident-wave  spectrum  and  the  weighting  factor  on  the 
incident  wave  spectrum  contribution  to  the  cost  function  is  set  at  <p  =  0.01.  The  iteration 
history  for  the  cost  function  J  is  shown  in  figure  3d)  where  it  is  seen  that  the  cost  function 
is  reduced  by  roughly  two  decades  in  ten  iterations  of  the  conjugate-gradient  procedure. 
Also  shown  in  figure  3d)  are  the  two  separate  terms  in  the  cost  function,  as  given  by  (28), 
where  J0bS  represents  the  first,  ‘observation’,  term  and  represents  the  second,  ‘boundary’, 
term.  The  contribution  to  the  cost  function  from  the  error  in  the  predicted  spectrum  at  the 
observation  location  (J0bs)  is  reduced  by  four  orders  of  magnitude,  and  is  still  dropping  after 
ten  iterations. 

Figures  3a)  and  3b)  show  the  estimates  for  the  incident-wave  and  observed  directional 
spectra  obtained  from  the  assimilation  procedure  The  significant  wave  height  for  the  esti¬ 
mated  incident-wave  spectrum  is  Hs  =  3.06  m.  about  3  percent  low,  and  the  peak  frequency 
fp  =  0.108  Hz  is  correct,  while  the  peak  wave  direction  is  9P  =  168.8°,  which  matches  the 
true  value  within  the  directional  resolution  of  the  computations.  The  under-estimation  of 
the  significant  wave  height  is  at  least  partly  due  to  the  introduction  of  the  incident  wave 
spectrum  E/,  into  the  cost  function;  this  serves  to  minimize  the  energy  in  the  incident-wave 
spectrum,  and  hence  puts  downward  pressure  on  the  wave  height.  Figure  3b)  shows  the 
estimated  directional  spectrum  at  the  FRF  8  m  array  location.  Here,  the  match  to  the 
spectrum  in  figure  2b)  is  almost  exact,  consistent  with  the  cost-function  contribution  J„b3 
in  figure  3d).  Also  shown  in  figure  3c)  is  the  estimated  significant- wave-height  field  which 
shows  excellent  agreement  with  the  result  shown  in  figure  2c). 

These  results  indicate  that  the  assimilation  procedure  is  capable  of  converging  to  a  solu¬ 
tion  that  is  almost  identical  to  those  from  the  forward  simulation  which  produced  the  data 
used.  The  results  shown  are  for  an  initial  guess  of  zero  for  the  incident  wave  spectrum.  An 
initial  guess  of  a  ‘white’  spectrum,  constant  for  all  wave  directions  and  frequencies,  has  also 
been  used.  This  produced  similar  results  to  the  zero-initial-guess  case,  and  demonstrates  the 
uniqueness  of  the  resulting  estimated  wave  spectra. 
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Figure  2.  Wave  directional  spectra  and  significant  wave  height:  a)  Incident-wave 
spectrum  applied  at  the  right-hand  (Eastern)  boundary  of  the  domain; 
b)  SWAN-calculated  spectrum  at  FRF  8  m  array  location;  c)  S WAN- 
calculated  significant  wave  height  Hs  for  region. 
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Figure  3.  Results  for  assimilation  of  synthetic  observations  at  the  FRF  8  m  array 
location:  a)  estimate  of  the  incident-wave  directional  spectrum;  b)  esti¬ 
mate  of  the  directional  spectrum  at  the  FRF  8  m  location;  c)  resulting 
significant  wave  height  for  region;  d)  iteration  history  for  the  cost  func¬ 
tion  J,  along  with  its  two  components  defined  in  (28). 
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Application  to  Wave-Spectrum  Observations 

The  foregoing  example  showed  that  the  assimilation  procedure  works  well  for  synthetic  data. 
To  demonstrate  its  potential  for  real-world  applications,  it  will  now  be  applied  to  a  measured 
directional-spectrum  data  set  from  the  FRF  8  m  array.  The  resulting  estimates  of  the  wave 
spectra  and  their  integral  properties  will  be  compared  to  both  those  of  the  8  m  array  spec¬ 
trum,  and  those  from  two  NDBC  measurement  stations  located  in  the  region:  NDBC  44014 
which  is  located  to  the  Northeast  of  the  observation  location  offshore  near  the  continental 
shelf  break,  and  CHLV2  near  the  entrance  to  the  Chesapeake  Bay,  in  shallower  water.  These 
NDBC  instruments  provide  frequency  spectra,  and  estimates  of  the  integral  parameters  of 
the  spectra. 

Figure  4  shows  the  results  for  assimilation  of  a  directional  spectrum  for  13  September  2001 
at  01:00  EST.  Here  again  (p  =  0.01  and  ten  conjugate-gradient  iterations  are  used.  For  this 
case  the  significant  wave  height  was  Hs  —  1.58  rn,  the  peak  frequency  was  fp  =  0.089  Hz 
and  the  mean  frequency  was  fm  =  0.130  Hz.  The  peak  wave  direction  was  0P  —  183.8°,  the 
mean  direction  at  the  peak  frequency  was  Bm,p  =  187.3°  and  the  mean  wave  direction  was 
()„,  1*8. IB  The  directional  spectrum  is  shown  in  figure  4a).  The  estimated  directional 
spectrum  from  the  issimilation  procedure  which  corresponds  to  this  observation  is  shown 
in  ti-up  lot  when  the  >ignifi<  ant  wave  height  is  lower  by  about  six  percent  (Hs  =  1.48  m) 
and  tit'1  wave  tire*  t i«>ii  bv  all  three  measures  used  above  is  estimated  within  two  degrees. 
The  peak  frequency  matches  the  observations,  but  the  mean  frequency  is  slightly  lower, 
indicating  thai  the  lower  wave  height  is  due  to  lower  energy  at  high  frequencies;  a  conclusion 
supported  by  comparison  of  figures  4a)  and  b).  Figure  4c)  shows  the  estimated  significant 
wave  height  field  for  the  region.  In  this  case,  the  waves  are  propagating  nearly  East  to  West 
and.  as  a  result,  there  are  triangular  regions  at  the  North  and  South  boundary  where  the 
wave-height  predictions  are  erroneously  low.  The  iteration  history  of  for  the  cost  function  J 
is  shown  in  figure  4d).  In  this  case  the  cost  function  decrease  is  slightly  more  than  an  order 
of  magnitude;  this  is  mainly  due  to  the  reduction  in  the  relatively  high  contribution  from 
the  observation  Jo(,s,  which  plateaus  after  about  seven  iterations.  This  is  indicative  of  the 
shortcomings  in  the  SWAN  model  physics  when  compared  to  real-world  data;  however,  the 
overall  agreement  as  shown  above  is  quite  good. 

A  more  stringent  test  of  the  assimilation  procedure  is  whether  the  estimated  wave  field 
matches  observational  data  other  than  that  used  as  input.  For  this,  contemporaneous  NDBC 
wave  data  is  used.  Figures  5a)  and  b)  show  the  direction  spectra  for  the  two  NDBC  locations, 
44014  and  CHLV2,  respectively. 

Figure  5c)  shows  a  comparison  of  the  wave  frequency  spectra  from  NDBC  44014  and 
SWAN,  where  the  SWAN  spectra  are  calculated  by  integrating  the  frequency-direction  spec¬ 
tra  in  6.  The  general  shape  of  the  spectrum  is  quite  well  captured  although  for  frequencies 
near  the  peak  and  lower,  SWAN  slightly  under-estimates  the  spectrum.  The  peak  frequency 
is  fp  =  0.090  Hz  from  the  data  while  SWAN  predicts  fp  =  0.081  Hz,  which  agrees  to  within 
the  resolution  of  the  two  discretizations  of  the  spectrum.  Comparing  integral  parameters,  it 
is  seen  that  SWAN  under-estimates  Hs  by  some  twenty  percent  ( Ha  =  1.59  m  vs.  2.00  m); 
this  is  due  to  the  under-estimation  of  the  spectrum  at  low  frequencies  and,  to  a  lesser  extent, 
to  the  fact  that  the  44014  spectrum  extends  to  higher  frequency  than  the  SWAN  spectrum 
grid.  The  mean  frequencies  agree  to  within  three  percent  (fm  =  1.37  Hz  for  44014,  and 
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Figure  4.  Results  for  assimilation  of  a  measured  wave  spectrum  at  the  FRF  8  m 
array  location  for  13  September  2001  at  01:00  EST:  a)  measured  FRF 
directional  spectrum;  b)  estimate  of  the  FRF  directional  spectrum;  c)  re¬ 
sulting  estimate  of  the  significant  wave  height  Hs  for  region;  d)  iteration 
history  for  the  cost  function  J,  along  with  its  two  components  defined 
in  (28). 
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1.34  Hz  for  SWAN).  The  mean  wave  direction  at  the  peak  frequency  is  9mtP  =  185.0°  for 
44014  and  184.5°  for  SWAN.  Overall  the  agreement  is  excellent  and,  since  NDBC  44014  is 
located  off-shore  near  the  shelf-break,  this  indicates  that  the  assimilation  procedure  does  a 
good  job  in  estimating  the  incident-wave  spectrum 

Figure  5d)  shows  a  similar  comparison  of  the  wave  frequency  spectra  from  NDBC  CHLV2 
and  SWAN.  Here,  the  shape  of  the  spectrum  is  reasonably  well-captured;  however,  the  peak 
is  over-estimated  by  a  factor  of  three.  This  is  balanced  by  under-estimates  of  the  energy 
density  at  higher  frequency  (and  the  extreme  low-frequency  tail)  to  produce  a  significant- 
wave-height  estimate  of  Hs  =  1.41  m  from  SWAN,  compared  to  1.40  m  for  CHLV2.  The  peak 
frequency  is  fp  =  0.100  Hz  from  CHLV2,  compared  to  the  SWAN  estimate  of  0.081  Hz.  The 
mean  frequency  fm  is  again  lower  for  SWAN  (0.120  Hz),  than  for  CHLV2  (0.169  Hz);  this 
is  due  to  the  observed  under-estimation  of  the  spectrum  at  high  frequencies.  Comparison  of 
wave  direction  estimates  are  not  possible,  since  CHLV2  only  provides  frequency  information. 

From  the  foregoing,  it  is  clear  that  the  comparisons  of  the  assimilation  results  to  data 
from  NDBC  buoy  44014  and  station  CHLV2  produce  mixed  results.  The  agreement  is  very 
good  for  44014,  while  for  CHLV2,  it  is  somewhat  less  satisfying.  Examination  of  the  locations 
of  44014  and  CHLV2  relative  to  the  FRF  8  m  array  and  CHLV2  (see  e.g.  figure  5c)  and  the 
fact  that  the  wave  direction  indicated  by  the  off-shore  buoy  44014  is  &m>p  =  185°.  would 
seem  to  indicate  that  if  the  spectrum  is  accurately  estimated  at  the  location  of  44014,  it 
would  be  reasonably  accurate  at  both  the  FRF  location  and  CHLV2.  Good  agreement  with 
the  data  is  found  at  the  FRF  location,  but  not  CHLV2.  The  reason  for  this  is  apparent 
in  figure  5c),  which  shows  a  drop  in  significant  wave  height  in  the  vicinity  of  CHLV2,  and 
to  the  North.  CHLV2  appear  to  be  located  in  the  region  near  the  upper  boundary  where 
the  wave-height  will  be  under-estimated  due  to  boundary  effects.  Be  that  as  it  may,  the 
key  to  obtaining  an  accurate  regional  estimate  of  the  wave  field,  is  accurately  determining 
the  off-shore,  incident-wave  spectrum.  The  excellent  results  obtained  in  the  comparison  to 
NDBC  44014,  indicates  the  promise  of  this  approach. 

4.2  Estimation  of  the  Wave  Field  from  SAR  Observations 

We  now  apply  the  assimilation  procedure  to  estimate  wave  spectra  from  SAR-image  data,  in 
particular,  to  SAR-image-spectrum  data.  We  will  focus  on  SAR-image  data  for  the  region 
around  the  USACE  Field  Research  Facility  in  Duck  NC.  Figure  6  shows  an  ERS  SAR  image 
overlaid  on  the  bathymetry  contours  for  the  region  containing  the  FRF  site.  The  image  is 
roughly  100  km  on  a  side  and  it  reaches  from  the  land  out  to  the  edge  of  continental  shelf. 
The  image  is  oriented  to  reflect  the  SAR  look  direction  (the  cross-track  direction)  of  —  78.1°T 
and,  due  to  the  orbit  geometry,  the  satellite  passes  over  this  region  at  10:44  EST. 

Waves  are  often  apparent  in  SAR  images,  as  seen  in  Figure  7,  which  shows  a  close-up  of 
the  near-shore  region  in  the  SAR  image  shown  in  6.  Waves  propagating  cross-shore,  with 
crests  oriented  parallel  to  the  shore,  are  clearly  visible.  Also  apparent  is  speckle  noise  due 
to  the  coherent  nature  of  the  SAR  imaging  process.  Speckle  is  a  major  source  of  noise  in 
the  SAR  assimilation  process.  Figure  8  shows  the  directional  (k)  spectrum  of  a  1.6  km- 
square  subset  of  the  the  SAR  image.  The  wave  propagation  direction  is  clearly  indicated 
by  the  spectral  peak;  however,  since  the  SAR  image  is  essentially  a  snapshot,  there  is  a 
180°  propagation-direction  ambiguity.  Also  evident  is  point-to-point  variability  expected  in 
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Figure  5.  Results  for  assimilation  of  a  measured  wave  spectrum  at  the  FRF  8  m 
array  location  for  13  September  2001  at  01:00  EST  and  comparison  of 
to  NDBC  buoy  data:  a)  estimate  of  directional  spectrum  for  location  of 
NDBC  44014;  b)  estimate  of  directional  spectrum  for  location  of  NDBC 
CHLV2;  c)  comparison  of  estimated  one-dimension  directional  spectrum 
to  data  from  NDBC  44014;  d)  comparison  of  estimated  one-dimension 
directional  spectrum  to  data  from  NDBC  CHLV2. 
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Figure  6.  Example  of  an  ERS  SAR  image  overlaid  on  the  bathymetry  contours 
for  the  region  around  the  US  ACE  FRF;  -f-  indicates  the  locations  of 
the  SAR  spectrum  observations  used  for  assimilation.  The  SAR  look 
direction  is  -78.1°T  (the  +y  axis  points  North)  and  the  satellite  track 
is  oriented  at  191. 9°T. 


19 


©1997  i;SA 

Figure  7.  Magnified  view  of  near-shore  waves  for  the  SAR  image  shown  above  in 
figure  6. 


a  spectral  estimate  based  on  a  finite-size  region. 

Application  to  Simulated  SAR-Spectrum  Data 

To  evaluate  the  assimilation  procedure,  independent,  of  complications  due  to  errors  in  both 
the  wave-field  and  the  SAR  image  modeling,  it  is  first  applied  to  a  set  of  simulated  observa¬ 
tions.  The  wave  spectrum  calculated  using  the  SWAN  model  for  41  locations  in  (indicated 
by  the  crosses  in  figure  6,  above)  was  used  along  in  the  nonlinear  SAR  model  to  generate  the 
observations.  The  results  are  shown  in  figure  9.  The  cost  function  history  indicates  that  the 
cost  function  decreases  by  about  two  orders  of  magnitude  in  about  twenty-five  iterations, 
and  so  there  is  a  good  match  between  the  predicted  SAR  spectrum  and  the  simulated  obser¬ 
vations.  Since  there  is  no  noise  in  added  to  the  data,  the  errors  are  indicative  of  information 
lost  in  going  from  the  wave  spectrum  to  the  SAR-image  spectrum,  as  a  result  of  imaging  ef¬ 
fects.  Also  shown  in  figure  9  is  a  comparison  between  the  estimated  and  actual  wave  spectra 
at  the  location  of  the  FRF  8  m  array.  The  peak  frequency  is  well-estimated  and  the  peak 
wave  direction  is  within  about  7°.  The  significant  wave  height  is  under-estimated  by  nearly 
20  percent,  mainly  due  to  high-frequency  waves  which  are  lost  in  the  SAR-imaging  process. 
It  should  be  noted  that  this  is  best-case  performance,  but  it  indicates  that  the  assimila¬ 
tion  procedure  is  convergent  and  capable  of  estimating  the  wave  spectrum  from  SAR-image 
spectrum  observations. 
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Figure  8.  Sample  SAR-image  spectrum  from  the  near-shore  region.  Here  the  spec¬ 
trum  has  been  rotated  and  the  SAR  look  direction  is  aligned  with  the 
ky  axis. 

Application  to  SAR-Spectrum  Observations 

We  now  apply  the  assimilation  procedure  to  SAR-image  spectrum  data  from  ERS  imagery. 
For  each  of  the  results  shown  below,  41  1.6  km  square  (128x128  pixel)  subsets  were  extracted 
from  the  SAR  image  at  selected  SWAN  grid  points.  These  subsets  were  Fourier  transformed 
and  the  and  the  image  power  spectra  were  calculated.  Starting  from  an  assumed  zero,  the 
incident  wave  spectrum  was  estimated  via.  the  assimilation  procedure.  The  final  estimate  of 
incident  spectrum  is  that  which  produces  the  best-fit  of  the  estimated  SAR  spectra  (using 
equation  11)  to  the  observed  spectra,  calculated  from  the  imagery.  For  these  purposes,  the 
best  fit  was  achieved  when  the  cost  function  ceased  to  decrease,  and  levelled  off.  In  the  results 
below,  the  resulting  wave  field  estimate  is  evaluated  by  comparing  to  contemporaneous  data 
from  the  FRF  8  m-array. 

The  first  image  is  from  05  March  1997  at  10:44  EST;  this  is  the  SAR  imagery  shown 
above  in  figures  6  and  7.  Results  for  the  assimilation  procedure  are  shown  in  figure  10.  The 
cost  function  is  reduced  by  more  than  half,  not  quite  as  good  as  the  results  obtained  with 
simulated  data.  The  estimated  significant-wave-height  field  is  also  shown,  where  a  wedge 
of  low  wave  height  is  seen  near  the  Northern  boundary;  this  is  due  to  the  Southwesterly 
wave  propagation  direction.  The  estimated  wave  spectrum  for  the  location  of  the  FRF  8  m 
array  is  also  shown  in  figure  10,  along  with  the  contemporaneous  in-situ  spectrum.  The 
estimated  peak  wave  direction  and  the  peak  frequency  compare  favorably  with  the  in-situ 
data,  with  the  peak  frequency  within  five  percent  and  the  peak  direction  within  8°;  however, 
the  significant  wave  height  is  under-estimated  by  more  than  30  percent.  This  is  due  to  the 
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Figure  9.  Results  for  assimilation  of  simulated  SAR  data  for  the  region  around  the 
USACE  FRF;  upper  left:  cost  function  iteration  history;  upper  right:  es¬ 
timated  significant  wave  height;  lower  right:  estimated  spectrum;  lower 
left,  actual  simulated  spectrum.  Locations  of  SAR  spectrum  observa¬ 
tions  used  for  assimilation  are  shown  in  6. 
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smaller  directional  spread  and  frequency  width  of  the  estimated  spectrum.  This  is  due  to 
the  loss  of  high-frequency  (short  wavelength)  waves  wave  information  in  the  SAR  imagery, 
as  well  as  information  at  high  azimuth  (along-track)  wavenumbers. 

A  sample  of  the  spectrum  results  from  the  assimilation  procedures  are  shown  in  figure  11. 
These  results  are  for  a  location  near  the  shore.  The  figure  shows  the  estimated  wave  spec¬ 
trum,  the  SAR  spectrum  estimate  calculated  from  that  wave  spectrum,  the  actual  SAR 
spectrum  observation,  and  difference  between  the  estimated  and  observed  SAR  spectrum. 
This  latter  quantity  provides  the  input  into  the  adjoint  model.  From  the  figure  it  is  clear  that 
the  estimated  SAR  spectrum  compares  well  with  the  observed  spectrum,  the  main  difference 
results  from  the  speckle  noise  present  in  the  observed  spectrum. 

Results  for  assimilation  of  an  ERS  SAR  image  from  09  April  1997  are  shown  in  figure  12. 
For  these  data,  the  cost  function  drops  only  a  few  percent  over  eight  iterations  and  the  the 
estimated  significant  wave  height  is  only  0.4  m  or  less.  Comparison  of  the  estimated  and 
actual  wave  spectra  for  the  FRF  8  rn  array  shows  that  the  actual  wave  height  at  this  location 
is  1.38  m  vs.  the  estimate  of  0.3  rn.  The  reason  for  this  discrepancy  can  be  found  in  the 
observed  FRF  spectrum  the  dominant  wave  direction  is  209. 9°T,  73°  from  the  SAR  look 
direction,  and  the  waves  are  relatively  high  frequency  (Tp  <  6  sec).  Hence,  the  wave  energy  is 
mainly  outside  the  azimuth  pass-band  of  the  SAR-imaging  process.  It  is  interesting  to  note 
that  the  small  amount  of  energy  near  the  SAR  pass-band  is  captured  by  the  assimilation 
procedure.  This  is  an  example  where  the  information  lost  in  the  SAR  imaging  process  cannot 
be  retrieved  by  the  assimilation  procedure. 

Figure  13  shows  the  results  for  another  SAR  image  from  03  July  1996.  For  this  case, 
the  waves  are  out  of  the  Southeast  with  a  significant  wave  height  of  less  than  a  meter.  The 
assimilation  procedure  converges,  when  the  cost  function  is  reduced  by  about  one-third.  For 
the  FRF  8  m-array  location,  the  estimated  wave  height  is  0.54  m,  compared  to  the  actual 
wave  height  of  0.68  m,  the  estimated  peak  frequency  fp  =  0.108  Hz  is  within  five  percent 
of  the  actual,  and  the  estimated  peak  wave  direction  9P  =  272°T  differs  by  about  12°  from 
the  actual.  Comparing  the  directional  spectrum  plots  shows  that  the  under-estimation  of 
the  significant  wave  height  is  again  due  to  the  loss  of  wave  energy  at  higher  frequencies  and 
larger  angles  relative  to  the  SAR  look  direction.  The  discrepancy  in  the  estimated  peak  wave 
direction  seems  less  pronounced  if  the  width  of  the  spectral  peak  is  taken  into  account 
the  data  show  a  broad,  slightly  bi-modal  peak,  where  the  estimate  shows  a  corresponding 
broad,  but  uni-modal,  peak  in  in  the  same  location. 

The  final  result  is  from  14  May  1997  and  has  waves  again  propagating  from  the  Southeast 
with  a  comparable  significant  wave  height,  about  half  a  meter.  The  results  for  assimilation 
of  this  image  are  shown  in  figure  14.  The  cost  function  in  this  case  is  reduced  by  about 
40  percent  at  convergence  and  the  agreement  is  excellent  for  wave  height,  frequency  and 
direction. 

4.3  Estimation  of  Bottom  Friction  from  SAR  Observations 

The  above  results  for  estimation  of  the  wave  field  from  SAR  observations  indicates  that  there 
is  sufficient  information  contained  in  the  SAR  imagery  to  get  a  reliable  estimate  of  the  waves. 
We  wish  to  extend  this  estimation  capability  to  the  the  bottom  friction  coefficient 
We  do  this  using  a  similar  approach  to  the  above,  but  now  we  include  adjustments  in  the 
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Figure  10.  Results  for  ERS  SAR  data  for  05  March  1997  10:44  EST  for  the  region 
around  the  USACE  FRF;  x  indicates  the  locations  of  the  SAR  spec¬ 
trum  observations  used  for  assimilation  and  -Vindicates  the  location  of 
the  FRF  8  m  array. 
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Figure  11 .  Sample  detailed  results  from  the  assimilation  procedure  for  x  =18.9  km 
and  y  =159  km:  the  estimated  wave  spectrum,  the  estimated  SAR- 
image  spectrum,  the  observed  SAR-image  spectrum,  and  the  difference 
between  the  estimated  and  the  observed  SAR-image  spectra.  Here,  red 
is  positive,  blue  is  negative,  zero  is  white,  and  the  SAR  look  direction 
is  aligned  with  the  ky  axis. 
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Figure  12.  Results  for  ERS  SAR  data  for  09  April  1997  10:44  EST  for  the  region 
around  the  USACE  FRF;  x  indicates  the  locations  of  the  SAR  spec¬ 
trum  observations  used  for  assimilation  and  -^indicates  the  location  of 
the  FRF  8  m  array.  For  this  case,  the  wave  energy  falls  substantially 
outside  the  SAR  pass-band. 
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Figure  13.  Results  for  ERS  SAR  data  for  03  July  1996  10:44  EST  for  the  region 
around  the  USAGE  FRF;  x  indicates  the  locations  of  the  SAR  spec¬ 
trum  observations  used  for  assimilation  and  -hindicates  the  location  of 
the  FRF  8  m  array. 
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Figure  14.  Results  for  ERS  SAR  data  for  14  May  1997  10:44  EST  for  the  region 
around  the  USACE  FRF;  x  indicates  the  locations  of  the  SAR  spec¬ 
trum  observations  used  for  assimilation  and  Aindicates  the  location  of 
the  FRF  8  m  array. 
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bottom  friction  coefficient  calculated  using  (42)  at  each  iteration.  This  was  applied  the  SAR 
imagery  from  14  May  1997  assuming  an  initial  friction  coefficient  of  Cb  —  0.0015,  the  default 
for  the  SWAN  model.  Figure  15  shows  results,  the  r.m.s.  bottom  velocity  (used  in  the  friction 
model)  and  the  friction  coefficient  C&(x)  estimated  using  the  assimilation  procedure.  The 
assimilation  procedure  indicates  that  the  the  friction  coefficient  should  increase  by  about 
30%  as  the  shore  is  approached  in  the  region  of  the  SAR  image  observations.  Additional 
results  are  shown  in  figure  16,  where  it  is  seen  the  cost-function  reduction  is  similar  to  that 
achieved  above  when  bottom  friction  was  not  considered,  but  the  off-shore  significant  wave 
height  is  much  larger  than  before.  Even  though  the  off-shore  significant  wave  height  is  larger 
than  previously  estimated,  the  agreement  between  the  estimated  and  actual  FRF  spectra  is 
similar  to  that  obtained  without  consideration  of  the  bottom  friction. 

These  results  indicate  that  the  estimation  of  the  bottom  friction  coefficient  from  SAR 
imagery  using  variational  assimilation  is  feasible,  but  it  must  be  applied  more  broadly  to 
determine  the  robustness  and  accuracy  of  the  procedure. 


5  Summary  and  Conclusions 

A  variational  assimilation  procedure  for  near-shore  waves  using  the  SWAN  model  has  been 
described.  The  assimilation  procedure  can  be  applied  to  wave  spectrum  observations  or 
SAR-image  spectra  obtained  from  satellite-based  sensors.  For  SAR  assimilation,  the  SWAN 
model  is  augmented  with  the  Hasselmann  &  Hasselmann  [10]  nonlinear  model  relating  the 
SAR-image  spectrum  to  the  local  wave  spectrum.  The  assimilation  procedure  is  used  to 
estimate  the  incident-wave  spectrum  at  the  boundary,  and  alternatively  the  bottom  friction 
coefficient,  which  produces  the  best  fit  between  the  model  predictions  and  the  data. 

Results  for  wave-spectrum  observations  are  shown  for  synthetic  observational  data,  as 
well  as  actual  data.  Application  to  synthetic  data  shows  that  the  assimilation  methodology 
can  be  used  to  accurately  determine  a  spatially  uniform  incident-wave  spectrum  based  on 
a  single  near-shore  spectrum  observation.  For  real-world  spectrum  data  from  the  FRF  8  m 
array,  the  results  yield  a  good  fit  to  the  observational  data,  and  compare  reasonably  well  to 
NDBC  Buoy  observations. 

Results  for  SAR-image-spectrum  observations  are  also  shown  for  synthetic  observational 
data  and  actual  data.  Application  to  synthetic  data  shows  that  the  assimilation  methodology 
can  be  used  to  accurately  estimate  the  wave  spectrum  from  an  array  of  simulated  SAR-image 
spectra;  the  main  error  is  the  loss  of  some  high-frequency  wave  energy.  When  applied  to  an 
array  of  SAR-image  spectra  obtained  from  ERS  satellite-based  SAR.  imagery,  the  results  are 
similarly  good,  but  when  the  waves  are  propagating  at  large  angles  relative  to  the  SAR  look 
direction,  the  wave  energy  does  not  appear  in  the  SAR-image  spectrum  due  to  well-known 
SAR-irnaging  effects,  and  cannot  be  retrieved. 

The  use  of  SAR  assimilation  for  the  estimation  of  the  bottom  friction  coefficient  and 
resulting  energy  dissipation  was  investigated  and  the  results  appear  to  be  promising,  but 
more  extensive  application  is  required  to  determine  the  robustness  and  overall  validity  of  the 
estimation  procedure. 
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Figure  15.  Results  for  estimation  of  bottom  friction  from  ERS  SAR  data  for  14 
May  1997  10:44  EST:  r.m.s.  bottom  orbital  velocity  used  in  the  fric¬ 
tion  model  and  estimated  bottom  friction  coefficient  C 5  showing  an 
increased  level  in  the  near-shore  region. 
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Figure  16.  Results  for  the  cost  function  and  wave  field  for  assimilation  of  ERS 
SAR  data  for  14  May  1997  10:44  EST  with  bottom-friction  estimation; 
x  indicates  the  locations  of  the  SAR  spectrum  observations  used  for 
assimilation  and  -(-indicates  the  location  of  the  FRF  8  m  array. 
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